Chainer チュートリアル 〜機械学習とデータ分析入門〜
7. 単回帰分析と重回帰分析
本章では、基礎的な機械学習手法として代表的な単回帰分析と重回帰分析の仕組みを、数式を用いて説明します。
よっしゃ、やってくぞ
本チュートリアルの主題であるディープラーニングの前に、単回帰分析と重回帰分析を紹介することには 2 つの理由があります。
理由述べてくれるのありがてぇ〜
1 つ目は、単回帰分析と重回帰分析の数学がニューラルネットワーク含めたディープラーニングの数学の基礎となるためです。
なるほど
2 つ目は、単回帰分析のアルゴリズムを通して微分、重回帰分析のアルゴリズムを通して線形代数に関する理解を深めることができるためです。
はい
機械学習手法は、教師あり学習 (supervised learning) 、教師なし学習 (unsupervised learning) 、強化学習 (reinforcement learning) に大別され、単回帰分析は教師あり学習に含まれます。 本チュートリアルで扱う多くの手法が教師あり学習です。
教師あり学習の中でも典型的な問題設定は 2 つに大別されます。 与えられた入力変数から、10 や 0.1 といった実数値を予測する 回帰 (regression) と、「赤ワイン」、「白ワイン」といったカテゴリを予測する分類 (classification) の 2 つです。
ふむ、こういうことね
教師あり学習
回帰 = 実数値を予測する
分類 = カテゴリを予測する
教師なし学習
強化学習
で、今回は「回帰」の話
単回帰 = 1 つの入力変数から 1 つの出力変数を予測
重回帰 = 複数の入力変数から 1 つの出力変数を予測
回帰分析を行うアルゴリズムでは、以下の 3 ステップを順番に考えていきます。
Step 1 : モデルを決める
Step 2 : 目的関数を決める
Step 3 : 最適なパラメータを求める
データの前処理 (preprocessing) をひとつ紹介します。 データの 中心化 (centering) は、データの平均値が 0 になるように全てのデータを平行移動する処理を言います。
中心化によるデータ前処理の利点の一つとして、調整すべきパラメータを削減できることが挙げられます。中心化処理を行うことで切片を考慮する必要がなくなるため、データ間の関係性を表現する直線の方程式を、 $ y_c=wx_c のように、簡潔に表現可能となります。
なるほど
単回帰分析については
モデルを一次関数で定義
前処理で中心化して切片を削除
目的関数を二乗和誤差として定式化
微分して最適なパラメータ(=一次関数の傾き)を求める
数値例を実際に計算してみる
という流れ
重回帰分析については
モデルをM次元空間における線形式として定義: $ y = w_1x_1 + w_2x_2 + \dots + w_Mx_M + b
$ x_0=1, w_0=bと置くことで$ bを吸収: $ y = w_1x_1 + w_2x_2 + \dots + w_Mx_M + w_0x_0 = \bm{w}^T\bm{x}
$ M+1次元のパラメータベクトル$ \bm{w}と$ M+1次元の変数ベクトル$ \bm{x}
$ N個の変数ベクトルがある場合に$ n番目の変数ベクトルを$ \bm{x_n}とし、その時の予測値を$ y_nすると
$ y_n = \bm{w}^T\bm{x_n} となる
ここで$ N個の予測値を持つベクトル$ \bm{y}を考えてみると
$ \bm{y} = \begin{bmatrix}y_1 \\ \vdots \\ y_N \end{bmatrix} = \begin{bmatrix}\bm{w}^T\bm{x_1} \\ \vdots \\ \bm{w}^T\bm{x_N} \end{bmatrix} = \begin{bmatrix}\bm{x_1}^T\bm{w} \\ \vdots \\ \bm{x_N}^T\bm{w} \end{bmatrix} = \begin{bmatrix}\bm{x_1}^T \\ \vdots \\ \bm{x_N}^T \end{bmatrix}\bm{w} = \begin{bmatrix}x_{10} & x_{11} & \dots & x_{1M} \\ \vdots & \vdots & \ddots &\vdots \\ x_{N0} & x_{N1} & \dots & x_{NM} \end{bmatrix}\bm{w} = \bm{X}\bm{w}
上式では $ \bm{A}^T\bm{B} = \bm{B}^T\bm{A} の法則を使っている
ここで、N×M 行列 X は、各行が各データを表しており、各列が各入力変数を表しています。 このような行列を、デザイン行列 (design matrix) と呼びます。
ヘェ〜
目的関数を二乗和誤差として定式化
$ L = (t_1−y_1)^2+(t_2−y_2)^2+⋯+(t_N−y_N)^2 = \begin{bmatrix}t_1-y_1 && \dots && t_N-y_N \end{bmatrix}\begin{bmatrix}t_1-y_1 \\ \vdots \\ t_N-y_N \end{bmatrix} = (\bm{t}-\bm{y})^T(\bm{t}-\bm{y}) \\ = (\bm{t}-\bm{Xw})^T(\bm{t}-\bm{Xw})
目的関数を微分してパラメータベクトル$ \bm{w}を求める
そのために目的関数$ Lを微分しやすい形に色々変形してる
8. NumPy 入門
本章の目標は、単回帰分析と重回帰分析の章で学んだ重回帰分析を行うアルゴリズムをNumPy を用いて実装することです。
おす、やってく。
これ、並行してRubyの Numo::NArray で実装していったら強くなれるんじゃないか？と思いついた
思いついただけでやるかどうかは別
慣習的に、numpy にはしばしば np という別名が与えられます。 コード中で頻繁に使用するモジュールには、短い別名をつけて定義することがよく行われます。
それでは、numpy を np という名前で import してみましょう。
はーい
NumPy では ndarray というクラスで多次元配列を表現します
ndarray オブジェクトは shape という属性 （attribute） を持っており、その多次元配列の形 （shape） が保存されています。
ndarray の形は、要素が整数のタプルで表され、要素数はその多次元配列の次元数 （dimensionality, number of dimensions） を表します。
形は、その多次元配列の各次元の大きさを順に並べた整数のタプルになっています。
shapeが (3,) の場合は、「1次元目の大きさが3で、2次元目以降が存在しない多次元配列」ということか。これはつまり「要素数が3の1次元ベクトル」だね
np.zeros, np.ones, np.full , np.eye 便便便利〜
e という 4×5 行列を表す多次元配列から、1 行 2 列目の値を取り出すには
e[0, 1] インデックス指定で取り出す
Python のリストと同様にスライス表記 （slicing） を用いて選択したい要素を範囲指定することができます
code:py
# 4 x 5 行列 e の真ん中の 2 x 3 = 6 個の値を取り出す
なるほどなぁ
ndarray のデータ型は、dtype という属性に保存されています。
一度あるデータ型で定義した配列のデータ型を別のものに変更するには、astype を用いて変換を行います。
3×3 行列 a と 3 次元ベクトル b という大きさのことなる配列を定義して、それらを足してみましょう。
code:py
a = np.array([
])
c = a + b
c
形が同じ行列同士の場合と同様に計算することができました。
これは NumPy が自動的にブロードキャスト（broadcast）と呼ばれる操作を行っているためです。
これなぁ。前勉強会に行った時にも話題に上がったけど、意図しない挙動につながるケースないのかな？便利ではあると思うんだけど...
なんやかんやでRubyでも書いてみた
code:rb
pry(main)> b = Numo::NArray.[](1,2,3)
pry(main)> c = a + b
行列同士の要素ごとの四則演算は、通常は行列の形が同じでなければ定義できません。
ですよね
しかし、前節の最後では 3×3 行列に 3 次元ベクトルを足す計算が実行できました。
これが要素ごとの計算と同じように実行できる理由は、NumPy が自動的に 3 次元ベクトル b を 3 つ並べてできる 3×3 行列を想定し、a と同じ形に揃える操作を暗黙に行っているからです。 この操作を、ブロードキャストと呼びます。
a と同じ形に揃える操作を暗黙に行っている
暗黙に行なわずにエラー吐いてくれた方が堅牢じゃない？
行いたい算術演算が、大きい方の配列の一部に対して繰り返し行われることで実現されるため、実際に小さい方の配列のデータをコピーして大きい配列をメモリ上に作成することは可能な限り避けられます。
同じ内容の巨大な配列を使うときにメモリ効率が良いとか？
とはいえ、堅牢にやりたい人ができないのは微妙じゃないかなぁ。なんか他に理由があるのかな？
いまのところ、ブロードキャストをデフォルトじゃなくてオプショナルにしてほしいなぁ、と思ってる
形の異なる配列同士の計算がブロードキャストによって可能になるためにはルールがあります。
それは、「2 つの配列の各次元が同じ大きさになっているか、どちらかが 1 であること」です。
ほう
サンプルコードを読んだけどあんまりしっくりは来なかった
行列の要素ごとの積は * を用いて計算できました。 一方、通常の行列同士の積（行列積）の計算は、* ではなく、別の方法で行います。
はい
方法は 2 種類あります。
1つは、np.dot() 関数を用いる方法です。 np.dot() は 2 つの引数をとり、それらの行列積を計算して返す関数です。
np.dot(A, B)
もう 1 つは、ndarray オブジェクトが持つ dot() メソッドを使う方法です。
A.dot(B)
多次元配列に含まれる値の平均・分散・標準偏差・最大値・最小値といった統計値を計算する方法
code:py
x = np.random.randint(0, 10, (8, 10))
x.mean()
x.var()
x.std()
x.max()
x.min()
x.mean(axis=1)
code:rb
pry(main)> x = Numo::Int32.new(8,10).rand(10)
=> Numo::Int32#shape=8,10 [2, 6, 0, 1, 1, 9, 0, 2, 6, 2, 4, 6, 4, 6, 1, 6, 8, 4, 8, 0, 2, 7, 0, 7, 1, 0, 0, 6, 9, 6, 9, 2, 9, 6, 5, 2, 3, 2, 1, 0, 8, 3, 0, 0, 0, 2, 3, 5, 0, 7, 6, 6, 1, 7, 1, 4, 4, 6, 9, 6, 8, 3, 6, 4, 2, 2, 8, 3, 3, 0, pry(main)> x.mean
Did you mean? median
Intだと小数関連の計算はできないっぽいなぁ。Floatにキャストする必要がある
code:rb
pry(main)> y = Numo::Float32.cast(x)
=> Numo::SFloat#shape=8,10 [2, 6, 0, 1, 1, 9, 0, 2, 6, 2, 4, 6, 4, 6, 1, 6, 8, 4, 8, 0, 2, 7, 0, 7, 1, 0, 0, 6, 9, 6, 9, 2, 9, 6, 5, 2, 3, 2, 1, 0, 8, 3, 0, 0, 0, 2, 3, 5, 0, 7, 6, 6, 1, 7, 1, 4, 4, 6, 9, 6, 8, 3, 6, 4, 2, 2, 8, 3, 3, 0, pry(main)> y.mean
=> 3.9749999046325684
pry(main)> y.var
=> 8.429746627807617
pry(main)> y.std
pry(main)> y.stddev
=> 2.90340256690979
pry(main)> y.max
=> 9.0
pry(main)> y.min
=> 0.0
pry(main)> y.mean(1)
pry(main)> y.mean(axis: 1)
標準偏差は stdメソッドではなくて stddevメソッド。
meanの引数はそのまま入れてもいいし、キーワード引数 axis: を使っても良い
NumPy を用いた重回帰分析
やっときたか！
code:py
# Xの定義
X = np.array([
])
# データ数（X.shape0) と同じ数だけ 1 が並んだ配列 ones = np.ones((X.shape0, 1)) # concatenate を使い、1 次元目に 1 を付け加える
X = np.concatenate((ones, X), axis=1)
# 先頭に 1 が付け加わったデザイン行列
X
code:rb
pry(main)> design_x = Numo::Int32[
pry(main)* ]
pry(main)> ones = Numo::Int32.ones(design_x.shape0, 1) concatenate はクラスメソッドを使ったけど、インスタンスメソッドもある。ので、
code:rb
design_x = ones.concatenate(design_x, axis:1) # この書き方でもOK
目標値
code:py
# t の定義
code:rb
パラメータ$ {\bf w} = ({\bf X}^{{\rm T}}{\bf X})^{\rm -1}{\bf X}^{\rm T}{\bf t}
code:py
# Step 1
xx = np.dot(X.T, X)
# Step 2
xx_inv = np.linalg.inv(xx)
# Step 3
xt = np.dot(X.T, t)
# Step 4
code:rb
pry(main)> xx = design_x.transpose.dot(design_x)
linalg を使うには別途 numo-linalg というgemを入れる必要があった
9. scikit-learn 入門
本章では、この scikit-learn というライブラリを用いて、データを使ってモデルを訓練し、評価するという一連の流れを解説し、Chainer を使ったディープラーニングの解説に入る前に、共通する重要な項目について説明します。
よっしゃ
機械学習の様々な手法を用いる際には、データを使ってモデルを訓練するまでに、以下の 5 つのステップがよく共通して現れます。
Step 1：データセットの準備
Step 2：モデルを決める
Step 3：目的関数を決める
Step 4：最適化手法を選択する
Step 5：モデルを訓練する
ふむ
Step 1：データセットの準備
米国ボストンの 506 の地域ごとの住環境の情報などと家賃の中央値の情報を収集して作られた Boston house prices dataset というデータセット
チャールズ川の川沿いかどうか (0 or 1) 必要なのか？笑
チャールズ川の近くは住みやすいのか住みにくいのか。知らない。
最後の MEDV 以外の 13 個の指標から、MEDV を予測する回帰問題に取り組んでみましょう。
MEDV 物件価格の中央値 ね。OK
訓練用データセットとテスト用データセットに分割しましょう。
はい
どのように分けるかには色々な方法がありますが、単純に全体の何割かを訓練用データセットとし、残りをテスト用データセットとする、といった分割を行う方法はホールドアウト法 (holdout method) と呼ばれます。
ほ〜。他にも方法あるのかぁ。
scikit-learn では、データセットから指定された割合（もしくは個数）のデータをランダムに抽出して訓練用データセットを作成し、残りをテスト用データセットとする処理を行う関数が提供されています。
便便便利
Step 2 ~ 4：モデル・目的関数・最適化手法を決める
scikit-learn で重回帰分析を行う場合は、LinearRegression クラスを使用します。 sklearn.linear_model 以下にある LinearRegression クラスを読み込んで、インスタンスを作成しましょう。
code:py
from sklearn.linear_model import LinearRegression
# モデルの定義
reg_model = LinearRegression()
上記のコードは、前述の 2 〜 4 までのステップを行います。
簡単かよ...最強かよ...
Step 5：モデルの訓練
次にモデルの訓練を行います。 scikit-learn に用意されている手法の多くは、共通して fit() というメソッドを持っており、 再利用可能なコードが書きやすくなっています。
code:py
# モデルの訓練
reg_model.fit(x_train, t_train)
モデルの訓練が完了しました。
いや...簡単かよ...
重回帰分析では、重み w とバイアス b の２つがパラメータでした。 求まった重み w の値は model.coef_ に、バイアス b の値は model.intercept_ に格納されています。
reg_model.coef_, reg_model.intercept_
モデルの訓練が完了したら、精度の検証を行います。
そうだね
LinearRegression クラスは score() メソッドを提供しており、入力値と目標値を与えると訓練済みのモデルを用いて計算した決定係数 (coefficient of determination) という指標を返します。
あったなぁそんなの
これは、使用するデータセットのサンプルサイズを $ N、$ n 個目の入力値に対する予測値を $ y_{n}、目標値を $ t_n、そしてそのデータセット内の全ての目標値の平均値を $ \bar{t} としたとき、
$ R^{2} = 1 - \dfrac{\sum_{n=1}^{N}\left( t_{n} - y_{n} \right)^{2}}{\sum_{n=1}^{N}\left( t_{n} - \bar{t} \right)^{2}}
決定係数の最大値は 1 であり、値が大きいほどモデルが与えられたデータに当てはまっていることを表します。
はい
code:py
# 精度の検証
reg_model.score(x_train, t_train) # => 0.7645451026942549
訓練済みモデルを用いて訓練用データセットで計算した決定係数は、およそ 0.765 でした。
これは訓練用データを用いた評価だね
訓練が終わったモデルに、新たな入力値を与えて、予測値を計算させるには、predict() メソッドを用います。 訓練済みのモデルを使ったこのような計算は、推論 (inference) と呼ばれることがあります。
ここからは検証用データを用いた話だね
code:py
22.6 という目標値に対して、およそ 24.94 という予測値が返ってきました。
簡単だなぁ
それでは、訓練済みモデルの性能を、テスト用データセットを使って決定係数を計算することで評価してみましょう。
code:py
reg_model.score(x_test, t_test) # => 0.6733825506400171
訓練用データセットを用いて算出した値（およそ 0.765）よりも、低い値がでてしまいました。
そうだねぇ
訓練時に用いたデータに対してはよく当てはまっていても、訓練時に用いなかったデータに対しては予測値と目標値の差異が大きくなってしまう現象を、過学習 (overfitting) と言います。
はい
各ステップの改善
Step 1 の改善：前処理
前処理 (preprocessing) とは、欠損値の補完、外れ値の除去、特徴量選択、正規化などの処理を訓練を開始する前にデータセットに対して行うことです。
手法やデータに合わせた前処理が必要となるため、適切な前処理を行うためには手法そのものについて理解している必要があるだけでなく、使用するデータの特性についてもよく調べておく必要があります。
そうだね
今回のデータは、入力値の値の範囲が CRIM, ZN, INDUS, ... といった指標ごとに大きく異なっています。 そこで、各入力変数ごとに平均が 0、分散が 1 となるように値をスケーリングする標準化 (standardization) をおこなってみましょう。
scikit-learn では sklearn.preprocessing というモジュール以下に StandardScaler というクラスが定義されています。 今回は、これを用いてデータセットに対し標準化を適用します。
code:py
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(x_train)
# 標準化
x_train_scaled = scaler.transform(x_train)
x_test_scaled = scaler.transform(x_test)
# モデルの訓練
reg_model = LinearRegression()
reg_model.fit(x_train_scaled, t_train)
# モデルの評価
reg_model.score(x_train_scaled, t_train) # => 0.7645451026942549
reg_model.score(x_test_scaled, t_test) # => 0.6733825506400195
結果は変わりませんでした。
あ、変わんないんだ。
なんで？
べき変換をする別の前処理を適用し、再度同じモデルの訓練を行ってみましょう。
code:py
from sklearn.preprocessing import PowerTransformer
scaler = PowerTransformer()
scaler.fit(x_train)
# 標準化
x_train_scaled = scaler.transform(x_train)
x_test_scaled = scaler.transform(x_test)
# モデルの訓練
reg_model = LinearRegression()
reg_model.fit(x_train_scaled, t_train)
# モデルの評価
reg_model.score(x_train_scaled, t_train) # => 0.7859862563286062
reg_model.score(x_test_scaled, t_test) # => 0.7002856551689581
結果が改善しました。
なんで？
データが正規分布に基づいている方が、モデルの精度があがるっぽい
scaler.fit(x_train) で RuntimeWarning: divide by zero encountered in log 出たよ。
StandardScalerとPowerTransformerでは0に対する処理が違うのかもなぁ
パイプライン化
前処理用の scaler と 重回帰分析を行う reg_model は、両方 fit() メソッドを持っていました。 scikit-learn には、パイプラインと呼ばれる一連の処理を統合できる機能があります。 これを用いて、これらの処理をまとめてみましょう。
ん？リファクタリング？
「各ステップの改善」はモデルの評価値を改善するための各ステップでのアプローチ、という文脈で受け取っていたけど、そういうことではないのか？
code:py
from sklearn.pipeline import Pipeline
# パイプラインの作成 (scaler -> svr)
pipeline = Pipeline([
('scaler', PowerTransformer()),
('reg', LinearRegression())
])
# scaler および reg を順番に使用
pipeline.fit(x_train, t_train)
# モデルの評価
pipeline.score(x_train, t_train) # => 0.7859862563286062
pipeline.score(x_test, t_test) # => 0.7002856551689581
パイプライン化を行うことで、x_train_scaled のような中間変数を作成することなく、同じ処理が行えるようになりました。 これによってコード量が減らせるだけでなく、評価を行う前にテスト用データセットに対しても訓練用データセットに対して行ったのと同様の前処理を行うことを忘れてしまうといったミスを防ぐことができます。
確かにこりゃ便利だ
今回はやらない
10. CuPy 入門
CuPy は NumPy と高い互換性を持つ数値計算ライブラリです。 NumPy で提供されている多くの関数を NVIDIA GPU (Graphics Processing Unit) で実行することで簡単に高速化できるように設計されています。
GPU (graphics processing unit) は 3D グラフィックスの描画や、画像処理を高速に計算できるように設計された演算装置です。一方、一般的な計算で使用される CPU (central processing unit) は、幅広く様々な処理で利用されることを想定して設計されています。
GPU は、条件分岐を多用するような複雑な計算には向かない一方、行列計算のようなシンプルな計算を大量に並列処理する必要がある場合は、CPU よりもはるかに高速な場合があります。
はい
CuPy の基本的な使い方は NumPy とほとんど同じです。 そのため NumPy の使い方を知っていれば、パッケージ名を numpy から cupy に置き換えるだけで、多くの関数が NumPy と同じ使い方で利用できます。
はい
Chainer では、デフォルトでは内部の計算に NumPy が使用され、GPU の使用を宣言した場合においては CuPy が使用されます。
はい 完
11. Pandas 入門
Pandas はデータ操作によく用いられるパッケージであり、CSV などの一般的なデータ形式で保存されたデータの読み込みや、条件を指定しての一部データの抽出など、機械学習手法で取り扱うデータを整理するのに便利です。
ここで、データの特徴をおおまかに調べるために便利な df.describe() メソッドを実行してみましょう。
これ便利だな〜
次に値に対する条件を指定してデータの選択を行う方法を紹介します。
便利だね〜
結果は、各要素が条件を満たすか、満たさないかを表す True、False が各要素の位置に格納されたデータフレームやシリーズとなります。 これをマスク (mask) と呼ぶことがあります。
12. Matplotlib 入門
グラフの描画を行う際は Matplotlib が便利です。 Colab では標準で Matplotlib を使ってプロットを行うと描画結果がノートブック上に表示されます。
いろんな図が作れて便利ね〜
統計図の作成を簡単に行えるように Matplotlib をベースに作られたライブラリに seaborn というものがあります。